﻿using System;
using System.ComponentModel;
using System.Diagnostics;
using System.Collections;
using System.Runtime.Serialization;
using System.Text;
using System.Security.Permissions;
using System.Runtime.InteropServices;

namespace Heaven.Mathematics
{
	/// <summary>
	/// Четырехмерная матрица
	/// </summary>
	public struct Matrix4
	{
		#region Публичные поля

        /// <summary>
        /// Элементы матрицы
        /// </summary>
		public double M11, M12, M13, M14,
                     M21, M22, M23, M24,
                     M31, M32, M33, M34,
                     M41, M42, M43, M44;

		#endregion

		#region Конструкторы

		/// <summary>
		/// Инициализация матрицы поэлементно
		/// </summary>
		public Matrix4(
			double M11, double M12, double M13, double M14,
			double M21, double M22, double M23, double M24,
			double M31, double M32, double M33, double M34,
			double M41, double M42, double M43, double M44
			)
		{
            this.M11 = M11; this.M12 = M12; this.M13 = M13; this.M14 = M14;
            this.M21 = M21; this.M22 = M22; this.M23 = M23; this.M24 = M24;
            this.M31 = M31; this.M32 = M32; this.M33 = M33; this.M34 = M34;
            this.M41 = M41; this.M42 = M42; this.M43 = M43; this.M44 = M44;
		}

		/// <summary>
		/// Инициализация матрицы поэлементно
		/// </summary>
		/// <param name="elements">Массив вещественных одинарной точности (построчный перебор)</param>
		public Matrix4(float[] elements)
		{
			Debug.Assert(elements != null);
			Debug.Assert(elements.Length >= 16);

			M11 = elements[0]; M12 = elements[1]; M13 = elements[2]; M14 = elements[3];
			M21 = elements[4]; M22 = elements[5]; M23 = elements[6]; M24 = elements[7];
			M31 = elements[8]; M32 = elements[9]; M33 = elements[10]; M34 = elements[11];
			M41 = elements[12]; M42 = elements[13]; M43 = elements[14]; M44 = elements[15];
		}

        /// <summary>
		/// Инициализация матрицы заданием столбцов
		/// </summary>
		/// <param name="column1">Вектор, первый столбик</param>
        /// <param name="column2">Вектор, второй столбик</param>
        /// <param name="column3">Вектор, третий столбик</param>
        /// <param name="column4">Вектор, четвертый столбик</param>
        public Matrix4(Vector4 column1, Vector4 column2, Vector4 column3, Vector4 column4)
		{
			M11 = column1.X; M12 = column2.X; M13 = column3.X; M14 = column4.X;
			M21 = column1.Y; M22 = column2.Y; M23 = column3.Y; M24 = column4.Y;
			M31 = column1.Z; M32 = column2.Z; M33 = column3.Z; M34 = column4.Z;
			M41 = column1.W; M42 = column2.W; M43 = column3.W; M44 = column4.W;
		}

		/// <summary>
		/// Конструктор копирования
		/// </summary>
		public Matrix4(Matrix4 m)
		{
			M11 = m.M11; M12 = m.M12; M13 = m.M13; M14 = m.M14;
			M21 = m.M21; M22 = m.M22; M23 = m.M23; M24 = m.M24;
			M31 = m.M31; M32 = m.M32; M33 = m.M33; M34 = m.M34;
			M41 = m.M41; M42 = m.M42; M43 = m.M43; M44 = m.M44;
		}

		#endregion

		#region Константы

		/// <summary>
		/// Четырехмерная нулевая матрица из вещественных одинарной точности
		/// </summary>
		public static readonly Matrix4 Zero = new Matrix4(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);

		/// <summary>
        /// Четырехмерная единичная матрица из вещественных одинарной точности
        /// </summary>
		public static readonly Matrix4 Identity = new Matrix4(
			1,0,0,0,
			0,1,0,0,
			0,0,1,0,
			0,0,0,1
			);

		#endregion

		#region Статические арифметические методы

		/// <summary>
		/// Сложение двух матриц
		/// </summary>
		/// <param name="a">Слагаемое матрица 4x4</param>
        /// <param name="b">Слагаемое матрица 4x4</param>
        /// <returns>Новый экземпляр матрицы - сумма</returns>
		public static Matrix4 Add(Matrix4 a, Matrix4 b)
		{
			return new Matrix4(
				a.M11 + b.M11, a.M12 + b.M12, a.M13 + b.M13, a.M14 + b.M14,
				a.M21 + b.M21, a.M22 + b.M22, a.M23 + b.M23, a.M24 + b.M24,
				a.M31 + b.M31, a.M32 + b.M32, a.M33 + b.M33, a.M34 + b.M34,
				a.M41 + b.M41, a.M42 + b.M42, a.M43 + b.M43, a.M44 + b.M44
				);
		}

		/// <summary>
		/// Сложение матрицы и скаляра
		/// </summary>
		/// <param name="a">Матрица 4х4</param>
		/// <param name="s">Скаляр</param>
		/// <returns>Новый экземпляр матрицы - сумма</returns>
		public static Matrix4 Add(Matrix4 a, double s)
		{
			return new Matrix4(
				a.M11 + s, a.M12 + s, a.M13 + s, a.M14 + s,
				a.M21 + s, a.M22 + s, a.M23 + s, a.M24 + s,
				a.M31 + s, a.M32 + s, a.M33 + s, a.M34 + s,
				a.M41 + s, a.M42 + s, a.M43 + s, a.M44 + s
				);
		}

		/// <summary>
		/// Разница двух матриц
		/// </summary>
		/// <param name="a">Уменьшаемая матрица 4х4</param>
		/// <param name="b">Вычитаемая матрица 4х4</param>
		/// <returns></returns>
		/// <remarks>result[x][y] = a[x][y] - b[x][y]</remarks>
		public static Matrix4 Subtract(Matrix4 a, Matrix4 b)
		{
			return new Matrix4(
				a.M11 - b.M11, a.M12 - b.M12, a.M13 - b.M13, a.M14 - b.M14,
				a.M21 - b.M21, a.M22 - b.M22, a.M23 - b.M23, a.M24 - b.M24,
				a.M31 - b.M31, a.M32 - b.M32, a.M33 - b.M33, a.M34 - b.M34,
				a.M41 - b.M41, a.M42 - b.M42, a.M43 - b.M43, a.M44 - b.M44
				);
		}

		/// <summary>
		/// Разница матрицы и скаляра
		/// </summary>
		/// <param name="a">Умешьшаемое - матрица 4х4</param>
		/// <param name="s">Вычитаемое - скаляр</param>
		/// <returns>Новый экземпляр матрицы 4х4, содержащий разницу</returns>
		public static Matrix4 Subtract(Matrix4 a, double s)
		{
			return new Matrix4(
				a.M11 - s, a.M12 - s, a.M13 - s, a.M14 - s,
				a.M21 - s, a.M22 - s, a.M23 - s, a.M24 - s,
				a.M31 - s, a.M32 - s, a.M33 - s, a.M34 - s,
				a.M41 - s, a.M42 - s, a.M43 - s, a.M44 - s
				);
		}

        /// <summary>
        /// Произведение матрицы и скаляра
        /// </summary>
        /// <param name="a">Матрица</param>
        /// <param name="s">Скаляр</param>
        /// <returns>Новый экземпляр матрицы 4х4, содержащий произведение</returns>
        public static Matrix4 Multiply(Matrix4 a, double s)
        {
            return new Matrix4(
                a.M11 * s, a.M12 * s, a.M13 * s, a.M14 * s,
                a.M21 * s, a.M22 * s, a.M23 * s, a.M24 * s,
                a.M31 * s, a.M32 * s, a.M33 * s, a.M34 * s,
                a.M41 * s, a.M42 * s, a.M43 * s, a.M44 * s
                );
        }

	    /// <summary>
		/// Произведение двух матриц
		/// </summary>
		/// <param name="a">Матрица 4х4</param>
        /// <param name="b">Матрица 4х4</param>
        /// <returns>Новый объект матрицы 4х4, содержащий произведение</returns>
		public static Matrix4 Multiply(Matrix4 a, Matrix4 b)
		{
			return new Matrix4(
				a.M11 * b.M11 + a.M12 * b.M21 + a.M13 * b.M31 + a.M14 * b.M41,
				a.M11 * b.M12 + a.M12 * b.M22 + a.M13 * b.M32 + a.M14 * b.M42,
				a.M11 * b.M13 + a.M12 * b.M23 + a.M13 * b.M33 + a.M14 * b.M43,
				a.M11 * b.M14 + a.M12 * b.M24 + a.M13 * b.M34 + a.M14 * b.M44,

				a.M21 * b.M11 + a.M22 * b.M21 + a.M23 * b.M31 + a.M24 * b.M41,
				a.M21 * b.M12 + a.M22 * b.M22 + a.M23 * b.M32 + a.M24 * b.M42,
				a.M21 * b.M13 + a.M22 * b.M23 + a.M23 * b.M33 + a.M24 * b.M43,
				a.M21 * b.M14 + a.M22 * b.M24 + a.M23 * b.M34 + a.M24 * b.M44,

				a.M31 * b.M11 + a.M32 * b.M21 + a.M33 * b.M31 + a.M34 * b.M41,
				a.M31 * b.M12 + a.M32 * b.M22 + a.M33 * b.M32 + a.M34 * b.M42, 
				a.M31 * b.M13 + a.M32 * b.M23 + a.M33 * b.M33 + a.M34 * b.M43,
				a.M31 * b.M14 + a.M32 * b.M24 + a.M33 * b.M34 + a.M34 * b.M44,

				a.M41 * b.M11 + a.M42 * b.M21 + a.M43 * b.M31 + a.M44 * b.M41,
				a.M41 * b.M12 + a.M42 * b.M22 + a.M43 * b.M32 + a.M44 * b.M42, 
				a.M41 * b.M13 + a.M42 * b.M23 + a.M43 * b.M33 + a.M44 * b.M43,
				a.M41 * b.M14 + a.M42 * b.M24 + a.M43 * b.M34 + a.M44 * b.M44
				);
		}

		/// <summary>
		/// Трансформация передаваемого вектора матрицей
		/// </summary>
		/// <param name="matrix">Матрица 4х4, задающая трансформацию</param>
		/// <param name="vector">Исходный вектор</param>
		/// <returns>Новый объект вектора, содержащий трансформированный исходный</returns>
		public static Vector4 Transform(Matrix4 matrix, Vector4 vector)
		{
			return new Vector4(
				(matrix.M11 * vector.X) + (matrix.M21 * vector.Y) + (matrix.M31 * vector.Z) + (matrix.M41 * vector.W),
				(matrix.M12 * vector.X) + (matrix.M22 * vector.Y) + (matrix.M32 * vector.Z) + (matrix.M42 * vector.W),
				(matrix.M13 * vector.X) + (matrix.M23 * vector.Y) + (matrix.M33 * vector.Z) + (matrix.M43 * vector.W),
				(matrix.M14 * vector.X) + (matrix.M24 * vector.Y) + (matrix.M34 * vector.Z) + (matrix.M44 * vector.W));
		}
        
		/// <summary>
		/// Транспонирование матрицы
		/// </summary>
		/// <param name="m">Исходная матрица 4х4</param>
		/// <returns>Новый объект матрицы 4х4, содержащая результат - транспонированную матрицу</returns>
		public static Matrix4 Transpose(Matrix4 m)
		{
			Matrix4 t = new Matrix4(m);
			t.Transpose();
			return t;
		}

		#endregion

		#region Перегруженные методы

		/// <summary>
		/// Получение хасп-ключа для этого экземпляра
		/// </summary>
		/// <returns>32-битное знаковое целое</returns>
		public override int GetHashCode()
		{
			return 
				M11.GetHashCode() ^ M12.GetHashCode() ^ M13.GetHashCode() ^ M14.GetHashCode() ^
				M21.GetHashCode() ^ M22.GetHashCode() ^ M23.GetHashCode() ^ M24.GetHashCode() ^
				M31.GetHashCode() ^ M32.GetHashCode() ^ M33.GetHashCode() ^ M34.GetHashCode() ^
				M41.GetHashCode() ^ M42.GetHashCode() ^ M43.GetHashCode() ^ M44.GetHashCode();
		}

		/// <summary>
        /// Проверка эквивалентности этого объекта с другим
		/// </summary>
		/// <param name="obj">Объект для сравнения</param>
		/// <returns>Истина если объекты эквивалентны</returns>
		public override bool Equals(object obj)
		{
			if (obj is Matrix4)
			{
				Matrix4 m = (Matrix4)obj;
				return 
					(M11 == m.M11) && (M12 == m.M12) && (M13 == m.M13) && (M14 == m.M14) &&
					(M21 == m.M21) && (M22 == m.M22) && (M23 == m.M23) && (M24 == m.M24) &&
					(M31 == m.M31) && (M32 == m.M32) && (M33 == m.M33) && (M34 == m.M34) &&
					(M41 == m.M41) && (M42 == m.M42) && (M43 == m.M43) && (M44 == m.M44);
			}
			return false;
		}

		/// <summary>
		/// Конвертация матрицы в строку
		/// </summary>
		/// <returns>Строка, представляющая матрицу</returns>
		public override string ToString()
		{
			StringBuilder s = new StringBuilder();
			s.Append(String.Format( "|{0}, {1}, {2}, {3}|\n", M11, M12, M13, M14));
			s.Append(String.Format( "|{0}, {1}, {2}, {3}|\n", M21, M22, M23, M24));
			s.Append(String.Format( "|{0}, {1}, {2}, {3}|\n", M31, M32, M33, M34));
			s.Append(String.Format( "|{0}, {1}, {2}, {3}|\n", M41, M42, M43, M44));

			return s.ToString();
		}

		#endregion

		#region Методы

        /// <summary>
        /// Трансформация передаваемого вектора матрицей (вектор переводится в однородные координаты)
        /// </summary>
        /// <param name="vector">Исходный вектор</param>
        /// <param name="matrix">Матрица, с использовнаием которой происходит трансформация</param>
        /// <returns>Новый объект вектора, содержащий трансформированный исходный</returns>
        public static Vector3 Transform(Matrix4 matrix, Vector3 vector)
        {
            double w = (matrix.M14 * vector.X) + (matrix.M24 * vector.Y) + (matrix.M34 * vector.Z) + matrix.M44;
            return new Vector3(
                ((matrix.M11 * vector.X) + (matrix.M21 * vector.Y) + (matrix.M31 * vector.Z) + matrix.M41) / w,
                ((matrix.M12 * vector.X) + (matrix.M22 * vector.Y) + (matrix.M32 * vector.Z) + matrix.M42) / w,
                ((matrix.M13 * vector.X) + (matrix.M23 * vector.Y) + (matrix.M33 * vector.Z) + matrix.M43) / w);
        }

		/// <summary>
		/// Подсчет детерминанта матрицы
		/// </summary>
		/// <returns>Скаляр, значение</returns>
		public double Determinant()
		{
			/*double det = 0.0f;
			for (int col = 0; col < 4; col++)
			{
				if ((col % 2) == 0)
					det += this[0, col] * Minor(0, col).Determinant();
				else
					det -= this[0, col] * Minor(0, col).Determinant();
			}
			return det;*/

            double d = 
                M14 * M23 * M32 * M41 - M13 * M24 * M32 * M41 - M14 * M22 * M33 * M41 + M12 * M24 * M33 * M41 +
                M13 * M22 * M34 * M41 - M12 * M23 * M34 * M41 - M14 * M23 * M31 * M42 + M13 * M24 * M31 * M42 +
                M14 * M21 * M33 * M42 - M11 * M24 * M33 * M42 - M13 * M21 * M34 * M42 + M11 * M23 * M34 * M42 +
                M14 * M22 * M31 * M43 - M12 * M24 * M31 * M43 - M14 * M21 * M32 * M43 + M11 * M24 * M32 * M43 +
                M12 * M21 * M34 * M43 - M11 * M22 * M34 * M43 - M13 * M22 * M31 * M44 + M12 * M23 * M31 * M44 +
                M13 * M21 * M32 * M44 - M11 * M23 * M32 * M44 - M12 * M21 * M33 * M44 + M11 * M22 * M33 * M44;

            return d;
		}

		/// <summary>
		/// Подсчет сопряженной матрицы
		/// </summary>
		/// <returns>Новый экземпляр матрицы 4х4, содержащий сопряженную матрицу</returns>
		public Matrix4 Adjoint() 
		{
            // FIXME: подозрение, что функция неверно работает - в топку ее!

			Matrix4 result = new Matrix4();
			for (int row = 0; row < 4; row++)
			{
				for (int col = 0; col < 4; col++) 
				{
					if (((col+row) % 2) == 0)
						result[row, col] = Minor(col, row).Determinant();
					else
						result[row, col] = -Minor(col, row).Determinant();
				}
			}

			return result;
		}

		/// <summary>
		/// Создание новой матрицы 3х3 из этой без указанной строки и столбца
		/// </summary>
		/// <param name="row">Строка подлежащая исключению</param>
		/// <param name="column">Столбец подлежащий исключению</param>
		/// <returns>Матрица 3х3, содержащая полученный минор</returns>
		public Matrix3 Minor(int row, int column) 
		{
			int r = 0;
			Matrix3 result = new Matrix3();
			for (int iRow = 0; iRow < 4; iRow++) 
			{
				int c = 0;
				if (iRow != row) 
				{
					for (int iCol = 0; iCol < 4; iCol++) 
					{
						if (iCol != column) 
						{
							result[r,c] = this[iRow, iCol];
							c++;
						}
					}
					r++;
				}
			}
			return result;
		}

		/// <summary>
		/// Подсчет средней линии - сумма диагональных элементов
		/// </summary>
		/// <returns>Возвращает сумму диагональных элементов</returns>
		public double Trace()
		{
			return M11 + M22 + M33 + M44;
		}

		/// <summary>
		/// Транспонирование матрицы
		/// </summary>
		public void Transpose()
		{
            Swap(ref M12, ref M21);
            Swap(ref M13, ref M31);
            Swap(ref M14, ref M41);
            Swap(ref M23, ref M32);
            Swap(ref M24, ref M42);
            Swap(ref M34, ref M43);
        }

        void Swap(ref double a, ref double b)
        {
            double c = a;
            a = b;
            b = c;
        }

        /// <summary>
        /// Инвесия матрицы
        /// </summary>
        /// <param name="matrix">Source matrix</param>
        /// <returns>Инвертированная матрица</returns>
        public static Matrix4 Inverse(Matrix4 matrix)
        {
            return matrix.Inverse();
        }

        /// <summary>
        /// Инвесия матрицы
        /// </summary>
        /// <returns>Инвертированная матрица</returns>
        Matrix4 Inverse()
        {
            //return Adjoint() * (1.0f / Determinant());

            Matrix4 inverse = Matrix4.Zero;
            
            inverse.M11 = M23 * M34 * M42 - M24 * M33 * M42 + M24 * M32 * M43 - M22 * M34 * M43 - M23 * M32 * M44 + M22 * M33 * M44;
            inverse.M12 = M14 * M33 * M42 - M13 * M34 * M42 - M14 * M32 * M43 + M12 * M34 * M43 + M13 * M32 * M44 - M12 * M33 * M44;
            inverse.M13 = M13 * M24 * M42 - M14 * M23 * M42 + M14 * M22 * M43 - M12 * M24 * M43 - M13 * M22 * M44 + M12 * M23 * M44;
            inverse.M14 = M14 * M23 * M32 - M13 * M24 * M32 - M14 * M22 * M33 + M12 * M24 * M33 + M13 * M22 * M34 - M12 * M23 * M34;

            inverse.M21 = M24 * M33 * M41 - M23 * M34 * M41 - M24 * M31 * M43 + M21 * M34 * M43 + M23 * M31 * M44 - M21 * M33 * M44;
            inverse.M22 = M13 * M34 * M41 - M14 * M33 * M41 + M14 * M31 * M43 - M11 * M34 * M43 - M13 * M31 * M44 + M11 * M33 * M44;
            inverse.M23 = M14 * M23 * M41 - M13 * M24 * M41 - M14 * M21 * M43 + M11 * M24 * M43 + M13 * M21 * M44 - M11 * M23 * M44;
            inverse.M24 = M13 * M24 * M31 - M14 * M23 * M31 + M14 * M21 * M33 - M11 * M24 * M33 - M13 * M21 * M34 + M11 * M23 * M34;

            inverse.M31 = M22 * M34 * M41 - M24 * M32 * M41 + M24 * M31 * M42 - M21 * M34 * M42 - M22 * M31 * M44 + M21 * M32 * M44;
            inverse.M32 = M14 * M32 * M41 - M12 * M34 * M41 - M14 * M31 * M42 + M11 * M34 * M42 + M12 * M31 * M44 - M11 * M32 * M44;
            inverse.M33 = M12 * M24 * M41 - M14 * M22 * M41 + M14 * M21 * M42 - M11 * M24 * M42 - M12 * M21 * M44 + M11 * M22 * M44;
            inverse.M34 = M14 * M22 * M31 - M12 * M24 * M31 - M14 * M21 * M32 + M11 * M24 * M32 + M12 * M21 * M34 - M11 * M22 * M34;

            inverse.M41 = M23 * M32 * M41 - M22 * M33 * M41 - M23 * M31 * M42 + M21 * M33 * M42 + M22 * M31 * M43 - M21 * M32 * M43;
            inverse.M42 = M12 * M33 * M41 - M13 * M32 * M41 + M13 * M31 * M42 - M11 * M33 * M42 - M12 * M31 * M43 + M11 * M32 * M43;
            inverse.M43 = M13 * M22 * M41 - M12 * M23 * M41 - M13 * M21 * M42 + M11 * M23 * M42 + M12 * M21 * M43 - M11 * M22 * M43;
            inverse.M44 = M12 * M23 * M31 - M13 * M22 * M31 + M13 * M21 * M32 - M11 * M23 * M32 - M12 * M21 * M33 + M11 * M22 * M33;
            
            return inverse * (1.0f / Determinant());
        }

		#endregion

		#region Операторы сравнения

		/// <summary>
		/// Проверка на эквивалентность двух матриц
		/// </summary>
		/// <param name="a">Матрица 4х4, левосторонний операнд</param>
        /// <param name="b">Матрица 4х4, правосторонний операнд</param>
        /// <returns>Истина, если матрицы равны, иначе ложь</returns>
		public static bool operator==(Matrix4 a, Matrix4 b)
		{
			return ValueType.Equals(a,b);		
		}

		/// <summary>
		/// Проверка на отличность двух матриц 4х4
		/// </summary>
        /// <param name="a">Матрица 4х4, левосторонний операнд</param>
        /// <param name="b">Матрица 4х4, правосторонний операнд</param>
        /// <returns>Истина, если матрицы различны, иначе ложь</returns>
        public static bool operator!=(Matrix4 a, Matrix4 b)
		{
			return !ValueType.Equals(a,b);
		}

		#endregion

		#region Бинарные операторы

		/// <summary>
		/// Сложение двух матриц
		/// </summary>
        /// <param name="a">Матрица 4х4, левосторонний операнд</param>
        /// <param name="b">Матрица 4х4, правосторонний операнд</param>
        /// <returns>Новый объект матрицы 4х4, содержащий сумму</returns>
		public static Matrix4 operator+(Matrix4 a, Matrix4 b)
		{
			return Matrix4.Add(a,b);
		}

		/// <summary>
		/// Сложение матрицы и скаляра
		/// </summary>
        /// <param name="a">Матрица 4х4, левосторонний операнд</param>
        /// <param name="s">Скаляр, правосторонний операнд</param>
        /// <returns>Новый объект матрицы 4х4, содержащий сумму</returns>
        public static Matrix4 operator+(Matrix4 a, double s)
		{
			return Matrix4.Add(a,s);
		}

		/// <summary>
        /// Сложение скаляра и матрицы
        /// </summary>
        /// <param name="s">Скаляр, левосторонний операнд</param>
        /// <param name="a">Матрица 4х4, правосторонний операнд</param>
        /// <returns>Новый объект матрицы 4х4, содержащий сумму</returns>
        public static Matrix4 operator+(double s, Matrix4 a)
		{
			return Matrix4.Add(a,s);
		}

		/// <summary>
		/// Разница двух матриц
		/// </summary>
        /// <param name="a">Матрица 4х4, левосторонний операнд</param>
        /// <param name="b">Матрица 4х4, правосторонний операнд</param>
        /// <returns>Новый объект матрицы 4х4, содержащий разницу</returns>
        public static Matrix4 operator-(Matrix4 a, Matrix4 b)
		{
			return Matrix4.Subtract(a,b);
		}

		/// <summary>
		/// Разница матрицы и скаляра
		/// </summary>
		/// <param name="a">Матрица 4х4, уменьшаемое</param>
		/// <param name="s">Скаляр, вычитаемое</param>
        /// <returns>Новый объект матрицы 4х4, содержащий разницу</returns>
        public static Matrix4 operator-(Matrix4 a, double s)
		{
			return Matrix4.Subtract(a,s);
		}

		/// <summary>
		/// Произведение двух матриц
		/// </summary>
        /// <param name="a">Матрица 4х4</param>
        /// <param name="b">Матрица 4х4</param>
        /// <returns>Новый объект матрицы 4х4, содержащий произведение</returns>
        public static Matrix4 operator*(Matrix4 a, Matrix4 b)
		{
			return Matrix4.Multiply(a,b);
		}

		/// <summary>
		/// Умножение матрицы на вектор
		/// </summary>
        /// <param name="matrix">Матрица 4х4</param>
        /// <param name="vector">Четырехмерный вектор</param>
		/// <returns>Трансформированный вектор</returns>
		public static Vector4 operator*(Matrix4 matrix, Vector4 vector)
		{
			return Matrix4.Transform(matrix, vector);
		}

        /// <summary>
        /// Умножение матрицы на вектор
        /// </summary>
        /// <param name="matrix">Матрица 4х4</param>
        /// <param name="vector">Треххмерной вектор</param>
        /// <returns>Трансформированный вектор</returns>
        public static Vector3 operator *(Matrix4 matrix, Vector3 vector)
        {
            return Matrix4.Transform(matrix, vector);
        }

        /// <summary>
        /// Умножение матрицы на вектор
        /// </summary>
        /// <param name="matrix">Матрица 4х4</param>
        /// <param name="vector">Треххмерной вектор</param>
        /// <returns>Трансформированный вектор</returns>
        public static Vector3 operator *(Vector3 vector, Matrix4 matrix)
        {
            return Matrix4.Transform(matrix, vector);
        }

        /// <summary>
        /// Умножение матрицы на скаляр
        /// </summary>
        /// <param name="matrix">Матрица 4х4</param>
        /// <param name="scalar">Скаляр</param>
        /// <returns>Матрица</returns>
        public static Matrix4 operator *(Matrix4 matrix, double scalar)
        {
            return Matrix4.Multiply(matrix, scalar);
        }

		#endregion

		#region Индексаторы

        /// <summary>
        /// Доступ по индексу
        /// </summary>
        public double this[int index]
        {
            get 
			{
                switch (index)
                {
                    case 0: return M11;
                    case 1: return M12;
                    case 2: return M13;
                    case 3: return M14;
                    case 4: return M21;
                    case 5: return M22;
                    case 6: return M23;
                    case 7: return M24;
                    case 8: return M31;
                    case 9: return M32;
                    case 10: return M33;
                    case 11: return M34;
                    case 12: return M41;
                    case 13: return M42;
                    case 14: return M43;
                    case 15: return M44;

                    default: 
                        throw new InvalidOperationException(
                            "Попытка считать из несуществующего элемента " + 
                            index.ToString() + " матрицы 4х4");
                }
			}
			set 
			{
                switch (index)
                {
                    case 0: M11 = value; break;
                    case 1: M12 = value; break;
                    case 2: M13 = value; break;
                    case 3: M14 = value; break;
                    case 4: M21 = value; break;
                    case 5: M22 = value; break;
                    case 6: M23 = value; break;
                    case 7: M24 = value; break;
                    case 8: M31 = value; break;
                    case 9: M32 = value; break;
                    case 10: M33 = value; break;
                    case 11: M34 = value; break;
                    case 12: M41 = value; break;
                    case 13: M42 = value; break;
                    case 14: M43 = value; break;
                    case 15: M44 = value; break;

                    default:
                        throw new InvalidOperationException(
                            "Попытка записать в несуществующий элемент " + 
                            index.ToString() + " матрицы 4х4");
                }
			}			
		}

        /// <summary>
		/// Доступ по индексу строка - столбец
		/// </summary>
		public double this[int row, int column]
		{
			get 
			{
				return this[ row*4 + column ];
			}
			set 
			{				
				this[ row*4 + column ] = value;
			}			
		}

		#endregion

        #region Компьютерная графика

        /// <summary>
        /// Создание матрицы, задающей левостороннюю ортогональную проекцию
        /// </summary>
        /// <param name="width">Длина видимого объема в ближней плоскости отсечения</param>
        /// <param name="height">Высота видимого объема в ближней плоскости отсечения</param>
        /// <param name="near">Расстояние до ближней плоскости отсечения (ось аппликат)</param>
        /// <param name="far">Расстояние до дальней плоскости отсечения (ось аппликат)</param>
        /// <returns>Левосторонняя ортогональная проекция</returns>
        public static Matrix4 Orthogonal(double width, double height, double near, double far)
        {
            if (far == Double.PositiveInfinity)
                return Orthogonal(width, height, near);

            return new Matrix4( 2 * near / width, 0, 0, 0,
                                0, 2 * near / height, 0, 0,
                                0, 0, far / (far - near), 1,
                                0, 0, near * far / (near - far), 0);
        }

        /// <summary>
        /// Создание матрицы, задающей левостороннюю ортогональную проекцию 
        /// с бесконечной дальней плоскостью отсечения
        /// </summary>
        /// <param name="width">Длина видимого объема в ближней плоскости отсечения</param>
        /// <param name="height">Высота видимого объема в ближней плоскости отсечения</param>
        /// <param name="near">Расстояние до ближней плоскости отсечения (ось аппликат)</param>
        /// <returns>Левосторонняя ортогональная проекция</returns>
        public static Matrix4 Orthogonal(double width, double height, double near)
        {
            return new Matrix4(
              2 * near / width, 0, 0, 0,
              0, 2 * near / height, 0, 0,
              0, 0, 1, -near,
              0, 0, 1, 0);
        }

        /// <summary>
        /// Создание матрицы, задающей правосторннюю? перспективную проекцию
        /// </summary>
        /// <param name="fovY">Угол зрения в вертикальной плоскости, в радианах</param>
        /// <param name="ratio">Отношение ширины окна картинной плоскости к его высоте</param>
        /// <param name="near">Расстояние до ближней плоскости отсечения (ось аппликат)</param>
        /// <param name="far">Расстояние до дальней плоскости отсечения (ось аппликат)</param>
        /// <returns>Правосторннаяя перспективная проекция</returns>
        public static Matrix4 Frustum(double fovY, double ratio, double near, double far)
        {
            if (far == float.PositiveInfinity)
                return Frustum(fovY, ratio, near);

            double h = 1.0f / (float)Math.Tan(fovY / 2.0f);
            double w = h / ratio;

            return new Matrix4(w, 0, 0, 0,
                               0, h, 0, 0,
                               0, 0, -far / (far - near), -1,
                               0, 0, -near * far / (far - near), 0);


        }

        /// <summary>
        /// Создание матрицы, задающей правостоннюю? перспективную проекцию, 
        /// где дальняя плоскость отсечения является бесконечной
        /// </summary>
        /// <param name="fovY">Угол зрения в ретикальной плоскости, в радианах</param>
        /// <param name="ratio">Отношение ширины окна картинной плоскости к его высоте</param>
        /// <param name="near">Расстояние до ближней плоскости отсечения (ось аппликат)</param>
        /// <returns>Правосторонняя перспективная проекция</returns>
        public static Matrix4 Frustum(double fovY, double ratio, double near)
        {

            double h = 1.0f / (float)Math.Tan(fovY / 2.0f);
            double w = h / ratio;

            return new Matrix4(
              w, 0, 0, 0,
              0, h, 0, 0,
              0, 0, 1, -near,
              0, 0, 1, 0);
        }

        /// <summary>
        /// Создание правостонней? матрицы вида, 
        /// задающей положение и направление глаз наблюдателя
        /// </summary>
        /// <param name="eye">Привязка вида (положение глаз наблюдателя)</param>
        /// <param name="at">Точка визирования (положение точки на которую ориентирована камера)</param>
        /// <param name="up">Вектор, направленый вверх (не обязательно {0.0f, 1.0f, 0.0f})</param>
        /// <returns>Правосторонняя матрица вида</returns>
        public static Matrix4 LookAt(Vector3 eye, Vector3 at, Vector3 up)
        {
            Vector3 zAxis = Vector3.Normalize(eye - at);
            Vector3 xAxis = Vector3.Normalize(Vector3.CrossProduct(up, zAxis));
            Vector3 yAxis = Vector3.Normalize(Vector3.CrossProduct(zAxis, xAxis));

            return new Matrix4(xAxis.X, yAxis.X, zAxis.X, 0,
                                xAxis.Y, yAxis.Y, zAxis.Y, 0,
                                xAxis.Z, yAxis.Z, zAxis.Z, 0,
                                -Vector3.DotProduct(xAxis, eye),
                                -Vector3.DotProduct(yAxis, eye),
                                -Vector3.DotProduct(zAxis, eye), 1);
        }

        /// <summary>
        /// Поворот вокруг оси X
        /// </summary>
        /// <param name="yaw">Угол в радианах</param>
        /// <returns>Правосторонняя матрица поворота</returns>
        public static Matrix4 RotateYaw(double yaw)
        {
            double cos = (float)Math.Cos(yaw);
            double sin = (float)Math.Sin(yaw);

            return new Matrix4(1.0f, 0.0f, 0.0f, 0.0f,
                               0.0f, cos,  sin,  0.0f,
                               0.0f, -sin, cos,  0.0f,
                               0.0f, 0.0f, 0.0f, 1.0f);
        }

        /// <summary>
        /// Поворот вокруг оси Y
        /// </summary>
        /// <param name="pitch">Угол в радианах</param>
        /// <returns>Правосторонняя матрица поворота</returns>
        public static Matrix4 RotatePitch(double pitch)
        {
            double cos = (float)Math.Cos(pitch);
            double sin = (float)Math.Sin(pitch);

            return new Matrix4(cos, 0.0f, -sin, 0.0f,
                               0.0f, 1.0f, 0.0f, 0.0f,
                               sin,  0.0f, cos, 0.0f,
                               0.0f, 0.0f, 0.0f, 1.0f);
        }

        /// <summary>
        /// Поворот вокруг оси Z
        /// </summary>
        /// <param name="role">Угол в радианах</param>
        /// <returns>Правосторонняя матрица поворота</returns>
        public static Matrix4 RotateRole(double role)
        {
            double cos = (float)Math.Cos(role);
            double sin = (float)Math.Sin(role);

            return new Matrix4(cos, sin, 0.0f, 0.0f,
                               -sin, cos, 0.0f, 0.0f,
                               0.0f, 0.0f, 1.0f, 0.0f,
                               0.0f, 0.0f, 0.0f, 1.0f);
        }

        /// <summary>
        /// Матрица вращения, задается углами Эйлера
        /// </summary>
        /// <param name="yaw">Вокруг оси абсцис: yaw</param>
        /// <param name="pitch">Вокруг оси ординат: pitch</param>
        /// <param name="role">Вокруг оси аппликат: role</param>
        /// <returns>Правосторонняя матрица поворота</returns>
        public static Matrix4 Rotate(double yaw, double pitch, double role)
        {
            // FIXME: оптимизировать
            return RotateYaw(yaw) * RotatePitch(pitch) * RotateRole(role);

            /*
            // FIXME: похоже правосторонний вариант?
            double sa = Trigonometry.Sin(yaw);
            double ca = Trigonometry.Cos(yaw);
            double sb = Trigonometry.Sin(pitch);
            double cb = Trigonometry.Cos(pitch);
            double sg = Trigonometry.Sin(role);
            double cg = Trigonometry.Cos(role);
            double sbsa = sb * sa;
            double sbca = sb * ca;

            Matrix4 m = Matrix4.Zero;

            m.M11 = cg * cb;
            m.M12 = sg * cb;
            m.M13 = -sb;
            m.M14 = 0;
            m.M21 = cg * sbsa - sg * ca;
            m.M22 = sg * sbsa + cg * ca;
            m.M23 = cb * sa;
            m.M24 = 0;
            m.M31 = cg * sbca + sg * sa;
            m.M32 = sg * sbca - cg * sa;
            m.M33 = cb * ca;
            m.M34 = 0;
            m.M41 = 0;
            m.M42 = 0;
            m.M43 = 0;
            m.M44 = 1.0f;

            return m;*/
        }

        /// <summary>
        ///	Матрица масштабирования
        ///		|Sx 0  0  0 |
        ///		| 0 Sy 0  0 |
        ///		| 0  0 Sz 0 |
        ///		| 0  0 0  1 |
        /// </summary>
        public static Matrix4 Scale(double x, double y, double z)
        {
            return new Matrix4(x, 0.0f, 0.0f, 0.0f,
                               0.0f, y, 0.0f, 0.0f,
                               0.0f, 0.0f, z, 0.0f,
                               0.0f, 0.0f, 0.0f, 1.0f);
        }

        /// <summary>
        ///	Матрица смещения
        ///		| 1  0  0  0 |
        ///		| 0  1  0  0 |
        ///		| 0  0  1  0 |
        ///		| Tx Ty Tz 1 |
        /// </summary>
        public static Matrix4 Translate(double x, double y, double z)
        {
            return new Matrix4(1.0f, 0.0f, 0.0f, 0.0f,
                               0.0f, 1.0f, 0.0f, 0.0f,
                               0.0f, 0.0f, 1.0f, 0.0f,
                               x, y, z, 1.0f);
        }

        /// <summary>
        ///	Элементы, отвечающие за сдвиг
        ///		| 1  0  0  0 |
        ///		| 0  1  0  0 |
        ///		| 0  0  1  0 |
        ///		| Tx Ty Tz 1 |
        /// </summary>
        public Vector3 Translation
        {
            get
            {
                return new Vector3(this.M41, this.M42, this.M43);
            }
            set
            {
                this.M41 = value.X;
                this.M42 = value.Y;
                this.M43 = value.Z;
            }
        }

        /// <summary>
        /// Точка, на которую направлена камера
        /// </summary>
        public Vector3 LookAtVector
        {
            get
            {
                return new Vector3(M31, M32, M33);
            }
            set
            {
                M31 = value.X;
                M32 = value.Y;
                M33 = value.Z;
            }
        }

        /// <summary>
        /// Вектор, обозначающий верх для данной камеры
        /// </summary>
        public Vector3 UpVector
        {
            get
            {
                return new Vector3(M21, M22, M23);
            }
            set
            {
                M21 = value.X;
                M22 = value.Y;
                M23 = value.Z;
            }
        }

        /// <summary>
        /// Вектор направленный вбок для данной камеры
        /// </summary>
        public Vector3 RightVector
        {
            get
            {
                return new Vector3(M11, M12, M13);
            }
            set
            {
                M11 = value.X;
                M12 = value.Y;
                M13 = value.Z;
            }
        }

        /// <summary>
        /// Get a transformation matrix which transforms 
        /// objects along to the given direction 
        /// (the object must have original direction along YAxis(!))
        /// </summary>
        /// <param name="direction">Direction</param>
        /// <returns>Transformation</returns>
        public static Matrix4 TransformAlongTo(Vector3 direction)
        {
            Matrix4 matrix = Matrix4.Identity;
            if ((!Vector3.ApproxEqual(direction, Vector3.YAxis, 0.001)) && (!Vector3.ApproxEqual(direction, -Vector3.YAxis, 0.001)))
            {
                Vector3 firstVector = Vector3.CrossProduct(Vector3.YAxis, direction);

                firstVector.Normalize();
                // Получаем перпендикуляр к полученному ранее перепендикуляру
                Vector3 secondVector = Vector3.Normalize(Vector3.CrossProduct(firstVector, direction));
                // Получаем матрицу поворота
                matrix.M11 = firstVector.X; matrix.M12 = firstVector.Y; matrix.M13 = firstVector.Z;
                matrix.M21 = direction.X; matrix.M22 = direction.Y; matrix.M23 = direction.Z;
                matrix.M31 = secondVector.X; matrix.M32 = secondVector.Y; matrix.M33 = secondVector.Z;
                //matrix.M11 = firstVector.X; matrix.M21 = firstVector.Y; matrix.M31 = firstVector.Z;
                //matrix.M12 = direction.X; matrix.M22 = direction.Y; matrix.M32 = direction.Z;
                //matrix.M13 = secondVector.X; matrix.M23 = secondVector.Y; matrix.M33 = secondVector.Z;
            }
            else if (Vector3.ApproxEqual(direction, -Vector3.YAxis, 0.001))
            {
                matrix = Matrix4.Scale(1, -1, 1);
            }
            return matrix;
        }

        /*/// <summary>
        ///	Элементы, отвечающие за масштабирование
        ///		|Sx 0  0  0 |
        ///		| 0 Sy 0  0 |
        ///		| 0  0 Sz 0 |
        ///		| 0  0 0  0 |
        /// </summary>
        public Vector3 Scale
        {
            get
            {
                return new Vector3(this.M11, this.M22, this.M33);
            }
            set
            {
                this.M11 = value.X;
                this.M22 = value.Y;
                this.M33 = value.Z;
            }
        }*/

        // TODO: добавить извлечение компонентов вращения и 
        // мастабирования, но указать на трудоемкость операции

        #endregion
     }
}
